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ABSTRACT 

We  consider  bearing  estimation  of  multiple  narrow-band 
plane  waves  impinging  on  an  array  of  sensors.  For  this 
problem,  bearing  estimation  algorithms  such  as  minimum 
variance  distortionless  response  (MVDR),  multiple  signal 
classification,  and  maximum  likelihood  generally  require  the 
array  covariance  matrix  as  sufficient  statistics.  Interestingly, 
the  rank  of  the  array  covariance  matrix  is  approximately  equal 
to  the  number  of  the  sources,  which  is  typically  much  smaller 
than  the  number  of  sensors  in  many  practical  scenarios.  In 
these  scenarios,  the  covariance  matrix  is  low-rank  and  can 
be  estimated  via  matrix  completion  from  only  a  small  subset 
of  its  entries.  We  propose  a  distributed  matrix  completion 
framework  to  drastically  reduce  the  inter-sensor  communica¬ 
tion  in  a  network  while  still  achieving  near-optimal  bearing 
estimation  accuracy.  Using  recent  results  in  noisy  matrix 
completion,  we  provide  sampling  bounds  and  show  how  the 
additive  noise  at  the  sensor  observations  affects  the  recon¬ 
struction  performance.  We  demonstrate  via  simulations  that 
our  approach  sports  desirable  tradeoffs  between  communica¬ 
tion  costs  and  bearing  estimation  accuracy. 

Index  Terms —  Array  signal  processing.  Covariance 
analysis.  Direction  of  arrival  estimation.  Matrix  completion 

1.  INTRODUCTION 

Many  algorithms  for  narrow-band  direction  finding  with  sen¬ 
sor  arrays,  parameter  estimation  of  multiple  sinusoidal  sig¬ 
nals  superimposed  with  noise,  and  discriminating  multiple 
overlapped  echoes  require  the  covariance  matrix  of  the  obser¬ 
vations  as  sufficient  statistics  for  estimation.  Representative 
members  include  minimum  variance  distortionless  response 
(MVDR),  multiple  signal  classification  (MUSIC),  and  max¬ 
imum  likelihood  (ML)  [1].  Although  these  algorithms  have 
rigorous  estimation  guarantees,  they  have  an  important  disad¬ 
vantage  when  the  sensors  are  untethered:  centralized  process- 

Email:  {aew2@rice.edu,  volkan. cevher@epfl.ch};  This  work  was  sup¬ 
ported  by  the  grants  NSF  CCF-043n50,  CCF-0728867,  CNS-0435425, 
and  CNS-0520280,  DARPA/ONR  N66001-08-  1-2065,  ONR  N00014-07- 
1-0936,  N00014-08-1-1067,  N00014-08-1-1112,  and  N00014-08-1-1066, 
AFOSR  FA9550-07- 1-0301,  ARO  MURl  W311NF-07-1-0185,  and  the 
Texas  Instruments  Leadership  University  Program. 


ing.  Hence,  wireless  communications  among  the  sensors  can 
create  many  bottlenecks  and  critical  points  of  failure. 

In  a  practical  system,  estimating  the  array  covariance  ma¬ 
trix  in  a  distributed  fashion  can  provide  the  ability  to  scale,  re¬ 
duce  adversarial  vulnerability,  decrease  communication  costs, 
and  share  processing  responsibilities  among  individual  sen¬ 
sors  as  the  number  of  sensors  increase.  Distributed  processing 
over  the  network  implies  that  the  sensors  only  communicate 
with  their  immediate  neighbors  on  the  communication  graph 
with  a  message  passing  scheme  that  uses  constant  bandwidth. 
In  this  case,  the  greatest  challenge  becomes  finding  a  message 
passing  scheme  and  an  associated  algorithm  with  provable  es¬ 
timation  guarantees. 

In  this  paper,  we  propose  propagating  local  correlation  es¬ 
timates  over  the  edges  of  a  communication  graph  to  partially 
fill  in  the  entries  of  the  sensor  array’s  covariance  matrix.  We 
pay  particular  attention  to  tree  structured  graphs  generated  by 
Prim’s  algorithm  and  Dijkstra’s  algorithm  [2].  Both  of  these 
algorithms  take  in  a  set  of  nodes  and  construct  a  tree  struc¬ 
tured  graph  that  minimizes  a  cost  function,  such  as  one  based 
on  communication  costs.  We  further  propose  extensions  to 
this  simple  communication  model  by  greedily  propagating 
additional  signal  information  through  the  tree  in  order  to  in¬ 
crease  the  number  of  observed  matrix  entries  when  needed. 

For  the  bearing  estimation  problem,  we  discuss  why  the 
array  covariance  matrix  is  low  rank  for  most  practical  sce¬ 
narios.  Then,  we  exploit  recent  results  in  matrix  completion 
theory  to  recover  the  full  covariance  matrix  from  the  subset 
of  observed  entries  [3,4].  We  provide  analytic  bounds  on  the 
covariance  estimation  error  as  a  function  of  the  number  of 
time  snapshots,  the  additive  noise  variance,  and  the  number 
of  filled  matrix  entries.  While  doing  so,  we  restate  the  matrix 
sampling  bounds  in  [4]  within  the  context  of  bearing  estima¬ 
tion.  We  subsequently  investigate  the  accuracy  of  the  bear¬ 
ing  estimates  based  on  the  recovered  covariance  matrix  using 
MVDR  and  compare  the  results  with  the  estimates  based  on 
the  fully  sampled  covariance  matrix.  We  show  via  simulations 
that  our  approach  sports  desirable  tradeoffs  between  commu¬ 
nication  costs  and  bearing  estimation  accuracy. 

The  paper  is  organized  as  follows.  In  Section  2  we  provide 
the  necessary  background  on  signal  processing  for  bearing  es¬ 
timation  and  matrix  completion.  In  Section  3  we  provide  de¬ 
tails  of  our  distributed  estimation  approach  to  determine  the 
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array  covariance  matrix  and  provide  performance  bounds.  In 
Section  4  we  provide  simulation  results  to  demonstrate  the 
validity  of  our  approach  with  both  synthetic  and  real-world 
data.  We  conclude  in  Section  5. 

2.  BACKGROUND 

2.1.  Narrow-band  bearing  estimation  problem 

The  problem  of  bearing  estimation  of  multiple  narrow-band 
plane  waves  with  sensor  arrays  can  be  reduced  to  the  follow¬ 
ing  model  [1,5]: 

y{tn)  =  A{e)x{tn) +w{tn),  u  =  1, . . . ,  N;  (1) 

where  y  G  (complex)  is  the  noisy  data  vector  of  the 

sensor  network,  x  G  is  the  vector  of  unknown  signal 

amplitudes,  and  w  G  is  an  additive  noise. 

In  (1),  A{6)  G  has  the  following  structure 

A{e)  =  [  a{9i)  . . .  a{0K)  ]  ,  (2) 

where  a{9k)  G  is  called  the  steering  vector,  0^  G  M 

is  the  bearing  of  the  fc-th  source,  and  6  —  \^  9i  ...  9  k  \- 

We  denote  X  G  as  the  source  matrix,  where  X  = 

[x{ti), . . .  ,x(tN)]-  The  array  data  matrix  Y  G  and 

the  noise  matrix  W  G  are  similarly  defined. 

The  objective  of  the  bearing  estimation  problem  is  to  de¬ 
termine  9  given  the  noisy  observations  in  (1).  Without  loss 
of  generality,  we  focus  on  2-D  bearing  estimation;  hence,  the 
p-th  entry  in  the  steering  vector  a{9k)  for  the  fc-th  source  has 
the  following  well-known  expression: 

[a(0fc)]p  =  exp  {-jujpp  cos  {tpp  -  9k)}  ,  (3) 

where  {pp,  ipp)  is  the  p-th  sensor  position  in  polar  coordinates 
(angle  is  measured  with  respect  to  the  horizontal  axis)  and  w 
is  the  known  narrow-band  frequency  of  the  sources. 

In  the  sequel,  we  make  the  following  assumptions: 

A1  (Sensor  Network):  P  ^  K.  The  number  of  sensors  P 
is  large  compared  to  the  number  of  sources  K.  Note  that 
P  >  K  is  a  necessary  condition  for  the  uniqueness  of  the 
bearing  estimates. 

A2  (Source  Signals):  E{[at(fi)]^}  =  0,  E  { [a;(fi)]J  = 
0  (fc  ^  1),  and  E{[x{ti)]}  [x{tj)]i^}  =  0  (i  ^  j).  Source  sig¬ 
nals  are  independent  from  each  other  and  are  uncorrelated  in 
time.  The  symbols  '  and  E  denote  the  Hermitian  and  expec¬ 
tation  operators,  respectively. 

A3  (Noise):  [w{ti)]p  ~  Af  The  noise  has  an  inde¬ 

pendent  and  identically  distributed  (iid)  Gaussian  distribution 
with  zero  mean  and  known  variance 

2.2.  Covariance  based  bearing  estimation  algorithms 

Many  algorithms  to  estimate  9  exist  in  the  literature  [1,5]. 
The  defining  characteristics  of  these  algorithms  is  their  re¬ 
liance  on  the  the  sample  covariance  matrix  R,  given  below. 


as  sufficient  statistics: 

R=j^YY'  =  Z  +  a^I,  (4) 

where  Z  =  A(9)X X' A' (9) / N  is  a  rank  K  positive  semi- 
definite  matrix.  Example  algorithms  include  minimum  vari¬ 
ance  distortionless  response  (MVDR),  multiple  signal  char¬ 
acterization  (MUSIC),  and  maximum  likelihood  (ML). 

To  determine  the  source  bearings,  the  MVDR  algorithm 
calculates  a  power  vs.  bearing  pattern  via 

Pmvdr(0;  R)  =  [a\9)R-^a(9)]  ,  (5) 

whose  local  maxima  indicates  the  locations.  In  contrast,  the 
MUSIC  algorithm  relies  on  a  partial  eigenvector  expansion 
of  R  (in  a  similar  fashion  to  (5))  corresponding  to  the  noise 
subspace  by  assuming  K.  The  ML  solution  minimizes  the 
mean  squared  data  error,  which  can  be  directly  written  as  a 
function  of  J?  [  1  ] . 

2.3.  Matrix  completion  from  incomplete  data 

Here,  we  review  the  recent  results  on  noisy  matrix  completion 
that  will  be  necessary  to  quantify  the  performance  of  our  ap¬ 
proach.  We  start  by  assuming  that  we  only  observe  a  partial 
subset  of  the  entries  of  a  matrix  Z  (P  x  P).  Clearly,  in  the 
general  case,  we  would  be  at  a  loss  to  provide  any  meaningful 
estimate  of  the  missing  entries.  However,  in  the  special  case 
where  Z  is  of  low  rank  and  the  observations  are  noiseless, 
M  =  O  }P  log^  P)  randomly  chosen  entries  are  sufficient  to 
recover  all  the  entries  of  Z  with  probability  at  least  1  —  P~^ 
via  the  following  convex  optimization  problem  [3] 

Z  =  argmin  ||Z'||*  s.t.  Pa(Z)  =  Po(Z'),  (6) 

where  ||  •  ||*  is  the  nuclear  norm,  H  is  the  set  of  observed 
elements  in  the  matrix  Z,  and  Pa(Z)  retrieves  the  elements 
of  H  from  Z. 

In  this  paper,  we  are  directly  interested  in  the  case  when 
the  observed  entries  have  errors, 

[U].,  =  [Z],,  +  [E],,,  (7) 

where  [Ej^^  is  iid  zero  mean  Gaussian  with  variance  t^. 
Lor  this  case,  Candes  and  Yaniv  [4]  empirically  show  that 
stable  recovery  is  possible  when  the  size  of  the  set  H  is  a 
constant  times  the  degrees  of  freedom.  Lor  a  P  x  P  corre¬ 
lation  matrix  Z  of  rank  K,  the  number  of  required  entries 
is  0{K  (2P  —  AT)).  The  matrix  completion  is  then  accom¬ 
plished  by  solving: 

Z  =  argmin  i I |Pn(Z'  -  U)||i  +  p||Z'||*,  (8) 

where  |j  •  ||f  is  the  Lrobenius  norm.  Moreover,  the  relaxation 
parameter  is  chosen  as  y  =  yj2rnJPT  so  that  ||Pn(Z  — 

U)||f  <  5  where  5  =  s/m  +  s/SrriT,  where  m  =  |H|. 
The  estimation  error  can  then  be  approximately  bounded  by 

||Z-Z|1  <47|f5;see[4]. 


3.  DISTRIBUTED  COVARIANCE  ESTIMATION  ON 
GRAPHS 

Although  the  bearing  estimation  algorithms  discussed  in  Sec¬ 
tion  2  have  rigorous  estimation  guarantees  [1],  they  have  an 
important  disadvantage:  they  require  centralized  processing 
where  a  fusion  center  collects  Y  in  order  to  calculate  R.  The 
centralized  nature  of  these  algorithms  result  in  increased  com¬ 
munication  costs  with  many  bottlenecks  and  critical  points  of 
failure.  In  contrast,  we  seek  distributed  algorithms  since  they 
provide  the  ability  to  scale,  reduce  vulnerability,  decrease 
communication,  and  share  processing  responsibilities  among 
individual  sensors  as  the  number  of  sensors  increase. 

Suppose  now  we  are  given  a  communication  graph  Q, 
where  the  edges  correspond  to  communication  links  and  the 
vertices  correspond  to  the  sensors  with  some  pre-determined 
hierarchical  ordering  of  the  sensors,  e.g.,  a  routing  scheme. 
Distributed  processing  over  Q  implies  that  the  sensors  only 
communicate  with  their  immediate  neighbors  on  the  graph 
that  share  an  edge  with  a  message  passing  scheme  that  use 
constant  bandwidth.  The  graph  structure  Q  could  be  con¬ 
structed  with  multiple  objectives  in  mind,  such  as  minimiz¬ 
ing  communication  costs,  satisfying  a  minimum  connectiv¬ 
ity  probability,  etc.  In  general,  the  number  of  edges  is  much 
smaller  than  for  practical  purposes. 

3.1.  Covariance  estimation  on  trees 

To  demonstrate  the  ideas  more  concretely,  we  focus  on  tree 
structured  graphs  T,  obtained  by,  say.  Prim’s  Algorithm  or 
Dijkstra’s  algorithm  [2].  We  use  a  bottom-up  communication 
hierarchy  where  the  children  nodes  communicate  to  their  par¬ 
ents  their  observations  as  well  as  the  inherited  matrix  entries 
from  their  own  children.  In  this  case,  a  parent  sensor  can  im¬ 
mediately  compute  the  cross-correlation  terms  for  the  matrix 
R  corresponding  to  itself  and  its  children  nodes,  along  with 
all  pairwise  relations  between  children.  This  scheme  contin¬ 
ues  recursively  until  the  root  node  is  reached. 

It  is  easy  to  see  that  the  number  of  observed  entries  in  a 
covariance  matrix  is  lower  bounded  by  3P  —  2,  corresponding 
to  a  chain  graph.  Under  random  deployment.  Prim  or  Dijkstra 
algorithms  provide  much  higher  connectivity  than  this  lower 
bound  (empirically  enough  entries  for  two  targets).  Higher 
target  numbers,  however,  may  require  a  larger  number  of  en¬ 
tries  due  to  the  increased  rank  of  the  covariance  matrix.  When 
this  occurs,  signal  information  can  be  passed  up  through  mul¬ 
tiple  levels  of  the  tree  instead  of  only  to  the  parent.  Control¬ 
ling  the  amount  of  information  sent  through  the  tree  allows  us 
to  create  a  tradeoff  between  the  communication  cost  and  the 
number  of  observed  entries  in  the  covariance  matrix.  Passing 
all  signal  information  in  the  tree  reveals  all  entries  of  the  co- 
variance  matrix  but  incurs  the  maximum  communication  cost. 


3.2.  Recovering  the  full  covariance  matrix 

As  a  result  of  the  tree  message  passing  scheme,  the  root  node 
observes  a  partial  set  of  entries  H  of  the  covariance  matrix 
R  where  H  is  determined  by  the  edges  of  the  communication 
tree  T.  The  number  of  observed  entries  is  defined  by  m  = 
|H|.  The  following  lemmas  characterize  the  observed  entries 
of  R  and  the  recovered  covariance  matrix  R.  Proofs  of  these 
lemmas  can  be  found  in  [6]. 

Lemma  I  (Diagonal  entries  of  R).  Define  U  =  R  —  I. 
[U]pp  (p  =  I, .  ■ .  ,P)  is  approximately  zero  mean  Gaussian 
distributed  with  variance  <  4N~^  [^]pp 

Lemma  2  (Off-diagonal  entries  of  R).  (p,l  =  1,  •  •  • ,  T*; 
p  1)  is  approximately  zero  mean  Gaussian  distributed  with 
variance  -P  -f  cr^. 

Lemma  3  (Frobenius  norm  of  Rn  {E)).  Defined  =  \\Vci  (E)  ||/r. 
Then,  w.h.p.,  S  <  2aN~^/'^  maxp  [R]pp  {jn  +  s/Sr^. 

Lemma  4  (Estimation  guarantee  on  R).  With  S  defined  as 
above  and  with  E  =  R  —  R  we  have  the  following  recon¬ 
struction  guarantee,  |1R|1f  <  ^\J 

In  our  approach,  the  root  node  solves  the  optimiza¬ 
tion  problem  in  (8)  and  reports  R  -h  a^I.  We  use  = 
max{Tpp,Tpi)  as  characterized  in  the  lemmas  1  and  2  when 
choosing  the  relaxation  parameter  p,. 

We  can  now  bound  the  performance  of  the  beamformer 
operating  on  the  reconstructed  correlation  matrix.  We  fur¬ 
ther  assume  that  ||i?~^R|jF  <  1-  This  condition  implies 
that  R  is  sufficiently  well  conditioned  a  posteriori  via  di¬ 
agonal  loading,  which  is  common  in  the  application  of  the 
MVDR  beamformer.  We  denote  R'  —  R  -f  7I  as  the  reg¬ 
ularized  covariance  matrix.  As  shown  in  [6],  7  satisfying 

7  >  ||R||f  <  -I-  25  is  sufficient  to  guaran¬ 

tee  ||R~^R||f  <  1-  Often  lower  values  of  7  can  be  safely 
used.  Defining  k{R!)  as  the  condition  number  of  R!  we  can 
show  that  k{R!)  ~  ■ 

With  the  above  background,  we  can  state  the  following 
lemma  using  some  results  from  matrix  perturbation  theory 

[7]: 

Lemma  5.  Defining  k(R')  as  the  condition  number  of  R! 
and  assuming  that  ||R^~^.E||f  <  1-  Then  we  have  the  fol¬ 
lowing  bound  on  beamformer  deviation: 

PMVDRiO;  R')  -  Pmvdr{0;  R  )  <  Pmvdr{0;  R'). 

l-||R  E/IIf 

Our  experiments  with  synthetic  and  real  data  show  that 
this  worst-case  bound  is  somewhat  pessimistic  and  that  the  ac¬ 
tual  deviation  from  the  ideal  beampattern  is  often  quite  small. 
We  conjecture  that  the  errors  introduced  in  the  matrix  com¬ 
pletion  process  has  an  underlying  structure  that  should  be  ex¬ 
ploited  to  obtain  tighter  bounds. 


4.  SIMULATION  RESULTS 
4.1.  Communication  Tradeoffs 

As  explained  in  Section  3,  using  a  tree  based  routing  ap¬ 
proaches  in  conjunction  with  matrix  completion  we  can  form 
a  tradeoff  between  communication  costs  and  the  number  of 
entries  observed  in  our  autocorrelation  matrix,  R.  This,  in 
turn,  affords  a  tradeoff  between  the  overall  communication 
cost  and  the  resulting  accuracy  of  our  estimator.  Here,  we 
present  simulation  results  to  demonstrate  this  tradeoff  for  var¬ 
ious  values  of  the  array  signal  to  noise  ratio.  Each  simulation 
trial  uses  a  random  network  of  P  =  30  sensors  which  we 
link  together  via  Dijkstra’s  algorithm.  The  desired  number  of 
observed  matrix  entries  determines  how  much  information  is 
passed  up  the  tree.  Observing  the  full  matrix  corresponds  to 
passing  all  signals  up  to  the  fusion  center,  with  a  communi¬ 
cation  cost  equal  to  the  full  cost  of  the  Dijkstra  tree.  Each 
simulation  point  is  averaged  over  1000  trials. 


(a)  (b) 


Fig.  I.  Communication  and  accuracy  tradeoffs. 


4.2.  Bearing  estimation  experiments 

We  present  bearing  estimation  results  for  both  synthetic  and 
field  data.  Eor  synthetic  data,  we  simulate  two  targets  moving 
on  a  circular  track  and  estimate  bearing  angles  at  several  time 
steps.  We  deploy  a  random  array  of  P  =  30  and  connect 
the  network  via  Prim’s  algorithm.  We  then  pass  inter-sensor 
messages  as  described  in  Section  3  and  observe  200  out  of  a 
possible  900  matrix  entries.  The  root  node  performs  matrix 
completion  to  obtain  R  which  then  serves  as  the  input  to  its 
beamformer.  We  use  R'  —  R  +  O.ll  when  beamforming. 
Our  results,  displayed  in  Eigure  2,  show  that  bearing  estimates 
with  our  approach  are  quite  close  to  the  full  data  case. 


Fig.  2.  Two-target  bearing  tracks. 


Next,  we  demonstrate  our  approach  on  acoustic  held  data 
corresponding  to  a  multi-vehicle  convoy  moving  along  an  el¬ 
liptical  track  observed  by  a  network  of  10  sensors.  The  sen¬ 
sors  are  distributed  in  a  circle  of  3  m  diameter  with  the  refer¬ 
ence  sensor  in  the  center  of  the  ring.  Because  both  Dijkstra’s 
and  Prim’s  algorithm  would  force  full  connectivity  for  this 
deployment,  we  enforce  sparsity  as  shown  in  3(a)  resulting 
in  observing  40  out  of  the  100  entries  of  the  correlation  ma¬ 
trix.  No  diagonal  loading  was  used  in  this  example.  Estimate 
results  are  displayed  in  Eigure  3(b). 


(a)  Connectivity  graph  (b)  Experimental  comparison 

Fig.  3.  Field  data  results. 

5.  CONCLUSION 

We  discussed  distributed  low -rank  covariance  estimation  over 
tree-based  graphs  via  matrix  completion.  Our  approach  ex¬ 
ploited  the  ability  to  stably  reconstruct  low  rank  matrices 
from  highly  incomplete  partial  observations.  We  quantihed 
how  additive  white  Gaussian  noise  in  the  sensor  observations 
directly  affect  the  matrix  completion  performance.  We  ob¬ 
served  empirically  that  the  noise  is  well  controlled  and  has 
very  little  effect  on  the  bearing  estimation  of  the  MVDR 
beamformer.  As  future  work,  we  will  investigate  other  graph 
structures,  such  as  the  Bethe  lattice  [8]  and  loopy  graphs. 
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